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ABSTRACT 

Micro  air  vehicles,  especially  those  that  use  flapping  wings  for  propulsion  and  control,  are 
characterized  by  complex,  time-dependent,  nonlinear  interactions  amongst  a  variety  of  physical 
disciplines.  These  interactions  must  be  carefully  controlled  to  attain  desired  levels  of  performance, 
but  the  number  and  range  of  parameters  that  describe  possible  micro  aircraft  is  extremely  large, 
poorly  explored,  and  inadequately  understood.  Aerodynamic,  structural,  kinematic,  and  mechanical 
parameters  are  relevant  and  may  be  strongly  coupled  with  each  other.  In  this  report,  we  describe  a 
multidisciplinary  optimization  methodology  that  can  be  used  to  isolate  combinations  of  parameters 
that  elicit  beneflcial  aeroelastic  interactions  from  flapping  wings.  This  methodology  relies  on 
disciplinary  models  that  are  tailored  to  the  physics  of  each  problem  considered  and  on  highly 
efficient  solution  strategies  suitable  for  design  optimization.  We  then  report  on  results  from  the 
application  of  this  methodology  to  three  problems  involving  flapping  wings  in  either  hover  or 
steady  forward  flight.  In  each  case,  we  identify  physical  mechanisms  that  can  be  used  to  achieve 
significant  gains  in  relevant  performance  metrics  such  as  force  production  and  power  consumption. 
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INTRODUCTION 


Micro  air  vehicles  (MAVs)  are  small  autonomous,  semi-autonomous,  or  remotely  piloted 
aircraft.  Given  their  particularly  small  size  and  their  ability  to  operate  in  close  spaces,  they  are 
envisioned  to  perform  a  variety  of  missions  that  are  unavailable  or  inaccessible  to  conventionally 
sized  aircraft.  Intelligence,  surveillance,  and  reconnaissance  missions  are  of  particular  interest. 
MAVs  are  also  contemplated  as  micro  munitions,  capable  of  delivering  low-collateral  weapons 
into  restricted  spaces.  Many  other  missions  are  possible,  as  MAVs  have  the  potential  to  match 
the  variety  of  performance  characteristics  exhibited  in  the  natural  world  by  flying  insects,  birds, 
and  mammals.  However,  an  inadequate  understanding  of  the  complex,  nonlinear,  and 
multidisciplinary  physics  that  underlie  flight  at  these  scales  complicates  realization  of  such 
capabilities. 

Engineers  from  many  organizations  around  the  globe  have  already  successfully  constructed 
and  flown  a  wide  variety  of  small-scale  robotic  aircraft.  These  vehicles  range  in  size  from  small 
birds  upwards,  use  fixed-wing,  rotary,  or  flapping-wing  layouts,  and  exhibit  varying  levels  of 
endurance,  controllability,  and  agility.  AeroVironment,  Inc.  and  the  Defense  Advanced  Research 
Projects  Agency  (DARPA),  for  example,  have  recently  demonstrated  a  remotely  piloted  two- 
wing  flapping  MAV  weighing  19  g  with  a  wingspan  of  16  cm\  Called  the  "Nano 
Hummingbird,"  the  concept  demonstrator  is  capable  of  controlled  flight  and  hover  and  has  a 
sufficient  level  of  maneuverability  to  establish  a  desired  orientation  within  a  horizontal  plane,  to 
transition  between  hover  and  forward  flight,  and  to  fly  through  a  normal-sized  doorway.  The 
designers,  like  others  before  them,  met  their  flight  test  milestones  primarily  through  extensive 
experimental  testing  and  iterative  improvement  of  flight  test  hardware.  While  this  approach  was 
used  successfully  in  the  DARPA  program,  it  is  limited  by  the  degree  of  physical  insight  that  the 
designers  can  bring  to  the  workbench.  Given  the  multidisciplinary  and  nonlinear  character  of  the 
physics,  such  an  approach  is  expected  to  be  quite  limited.  We  believe  that  a  thorough 
understanding  of  aeroelastic  behavior  at  small  scales  and  an  associated  physics-based  design 
process  are  needed  to  realize  effective  MAV  designs. 

Engineers  engaged  with  a  MAV  design  will  be  driven  to  choose  a  structural  layout  and 
propulsive  mechanism  by  the  particulars  of  the  mission  requirements.  A  fixed  wing  design  will 
be  best  for  some  missions,  a  rotary  wing  design  for  other  missions,  and  a  flapping-wing  design 
for  still  others.  Quantitative  assessment  of  the  relative  strengths  of  each  configuration  type  is 
needed  to  help  guide  future  designs.  However,  we  have,  for  the  purposes  of  this  study,  chosen  to 
focus  on  flapping-wing  designs.  The  physics  of  flapping  flight  are  arguably  least  well 
understood,  and  flapping  wings  offer  significant  potential  for  agility. 

The  physics  of  flapping- wing  MAVs  is  challenging  due  to  the  nonlinear  and 
multidisciplinary  character  of  the  system.  As  is  well  known,  small  changes  in  a  system  parameter 
can  sometimes  produce  large  changes  in  the  system  response  of  a  nonlinear  system.  Also, 
couplings  between  two  or  more  physical  disciplines  can  give  rise  to  energy  transfers  that 
significantly  influence  system  outputs.  These  effects  are  often  deleterious  to  performance 
objectives,  but  can  occasionally  be  beneficial.  The  challenge  is  to  identify  those  parameters  that 
have  the  greatest  impact  on  the  system  performance,  exploit  beneficial  nonlinearities  and 
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multidisciplinary  interactions,  and  avoid  substantial  detrimental  effeets.  Multidisciplinary 
optimization  techniques  provide  a  systematie  approach  for  exploring  the  parameter  space  in 
systems  where  suffieiently  accurate  numerieal  models  are  available  to  represent  system  behavior. 
For  MAVs,  this  requires  an  aeroelastie  modeling  eapability  informed  by  both  theory  and 
experiment  that  can  capture  key  physieal  eouplings  between  the  unsteady  aerodynamies, 
structural  mechanics,  kinematics,  and  actuation  mechanics. 

There  are  many  ehallenges  assoeiated  with  modeling  flow  fields  around  small,  agile  aircraft. 
The  flow  physies  are  highly  nonlinear  and  inherently  dynamic.  In  fact,  MAVs  depend  upon 
unsteady  flow  meehanisms  to  generate  forces  sufficient  for  flight.  Operating  at  low  Reynolds 
numbers  of  10^  or  less,  viscosity  and  turbulence  can  be  important  faetors,  and  vortieal  flow 
features  are  dominant.  Wing  interaetions  with  the  fuselage  or  other  wings  can  be  strong  due  to 
the  vehicle  geometry  and  the  flapping  kinematies.  Aeroelastie  behavior  can  only  be  properly 
characterized  when  the  aerodynamics  are  modeled  aeeurately.  On  the  other  hand,  the  strength  of 
aeroelastie  interactions  depends  upon  a  large  number  of  system  eonfiguration  parameters 
defining  geometry,  structure,  kinematies,  controls,  and  other  eharacteristies.  Exploring  this  large 
parameter  space  is  eostly  and  requires  fast  simulations.  This  gives  rise  to  a  situation  in  which 
different  levels  of  aerodynamic  modeling  fidelity  are  appropriate.  Low  and  medium  fidelities  are 
needed  to  explore  the  parameter  space  and  identify  regions  of  interesting  physies.  High-fidelity 
methods  are  needed  to  fully  charaeterize  the  aerodynamies  associated  with  benefieial  physical 
interactions. 

Small-scale  aircraft,  particularly  flapping-wing  MAVs,  generally  use  very  thin  wings.  These 
might  be  eonstructed  of  plies  of  eomposite  material  or  of  an  even  thinner  flexible  membrane 
material.  Beams  or  battens  are  often  ineorporated  into  the  wings  to  provide  additional  stiffness. 
Struetural  deformations  in  sueh  wings  can  be  large,  especially  for  flapping,  and  must  be 
considered  to  be  nonlinear.  In  some  cases,  the  true  structure  can  be  reduced  to  a  set  of  equivalent 
nonlinear  beams.  In  general,  however,  MAV  wings  should  be  treated  as  nonlinear  plates  or 
membranes,  possibly  eombined  with  nonlinear  beam  elements  for  stiffening. 
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METHODOLOGY 


We  use  multidisciplinary  optimization  to  explore  flapping-wing  MAV  physics,  with  an  eye 
towards  leveraging  beneficial  interactions  that  may  exist  between  elements  of  the  flapping 
system.  In  particular,  we  have  been  concerned  with  multidisciplinary  couplings  between  the 
unsteady  aerodynamics,  nonlinear  structure,  flapping-wing  kinematics,  and  flapping  actuation 
mechanism.  Computational  optimization  requires  both  accurate  and  timely  models  for  each  of 
these  physical  disciplines.  Key  physical  interactions  and  the  sensitivity  of  those  interactions  to 
variations  in  select  system  parameters  must  be  captured  adequately  so  that  the  optimization 
process  can  reveal  useful  parameter  combinations.  This  drives  the  researcher  towards  high- 
fidelity  models.  On  the  other  hand,  when  considering  many  system  parameters,  the  optimization 
process  requires  a  large  number  of  potentially  costly  analyses  and  exploration  of  the  parameter 
space  can  become  intractable.  The  multidisciplinary  optimization  approach,  therefore,  requires 
that  the  physical  models  be  carefully  chosen  so  as  to  capture  the  necessary  physical  features 
while  minimizing  computational  cost.  These  decisions  are  problem  dependent,  so  it  is  useful  to 
have,  for  each  discipline,  a  collection  of  models  of  varying  fidelity  to  draw  upon.  This  allows  the 
investigative  method  to  be  tailored  to  the  particular  problem  under  consideration. 

This  section  will  describe  the  major  physical  models  and  numerical  methods  that  have  been 
developed  and/or  used  in  our  research.  The  following  section  will  show  applications  of  these 
various  methods  to  problems  relevant  to  flapping-wing  MAVs. 

A.  Aerodynamics  Methods 

The  unsteady  and  nonlinear  aerodynamics  of  flapping  wings  represent  perhaps  the  most 
challenging  aspect  of  aeroelastic  modeling  for  multidisciplinary  optimization.  Complicated  flow 
features  are  usually  present  and,  in  many  cases,  improving  performance  depends  on  wisely 
controlling  these  features.  This  suggests  that  high-fidelity  models  such  as  the  Reynolds-averaged 
Navier-Stokes  (RANS)  and  large  eddy  simulation  Navier-Stokes  (LES)  methods  are  best  suited 
to  MAV  applications.  However,  it  has  been  shown  through  comparisons  with  experiments  that, 
for  certain  problems,  lower-fidelity  aerodynamics  methods  are  able  to  capture  the  major  vortical 
flow  features  associated  with  flapping  while  requiring  significantly  reduced  computational 
resources.  These  lower-order  methods,  therefore,  have  been  used  extensively  in  multidisciplinary 
studies  of  flapping  physics. 

1.  Quasi-Steady  Blade  Element  Aerodynamics 

We  use  an  internally  developed  implementation  of  the  aerodynamic  model  discussed  by 
Berman  and  Wang^'"*.  Following  the  blade  element  assumption,  3D  lifting  surfaces  are  divided 
along  the  span  into  a  number  of  chordwise  segments.  The  total  force  is  obtained  by  summing  the 
forces  calculated  on  each  segment.  The  force  acting  on  an  individual  element  is  found  by  a 
model  formulated  to  study  free-falling  plates  and  is  equal  to  the  sum  of  four  components'’^. 


^  ^circulation  d"  ^added  mass  d"  ^viscous  dissipation  d"  ^rotating  frame 


(1) 
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The  expression  for  the  cireulation  and  viseous  terms  involve  translational  and  rotational  foree 
coeffieients  and  a  drag  coeffieient  that  are  typieally  taken  from  experimental  measurements  of 
wind  tunnel  models,  birds,  or  insects. 


X-force  (side)  Y-force  (thrust)  Z-force  (lift) 


Figure  1.  Quasi-steady  blade  element  model  verification:  force  generation  of  a  rigid 

flapping  wing. 


It  only  takes  seconds  to  execute  a  quasi-steady  blade  element  simulation,  yet  results  can  be  quite 
accurate  for  many  problems.  For  example.  Figure  1  shows  a  good  agreement  between  force 
generation  curves  obtained  by  the  blade  element  and  Navier-Stokes  methods  for  a  flapping  half¬ 
ellipse  planform.  The  wing,  with  a  50  mm  span  and  an  18.2  mm  mean  chord,  was  flapped  at  a 
frequency  of  26  Hz  with  the  following  prescribed  kinematics: 

TT 

4>  =  —cos(2nft) 

0  =  0  (2) 

n  ^  ^  n 

^  =  -^sinCZTT/t)  -h- 


2.  UVLM Aerodynamics 

The  unsteady  vortex  lattice  method  (UVLM)  offers  a  balance  between  the  computational 
efficiency  of  a  quasi-steady  blade  element  method  and  the  physical  accuracy  of  a  Navier-Stokes 
method.  The  equations  are  simpler  and  faster  to  solve  than  the  Navier-Stokes  equations.  While 
viscosity  is  ignored  and  flow  separation  must  be  modeled,  UVLM  is  able  to  capture  the  major 
vortical  flow  features  associated  with  dynamic  MAV  flows.  We  have  created  both  2D  and  3D 
UVLM  solvers^’^.  The  UVLM  is  well  suited  to  incorporating  models  for  turbulence,  and  we  have 
completed  this  integration  in  2D^. 

Karamcheti  and  Katz  and  Plotkin  described  the  basics  of  this  method*®’^*.  It  is  developed 
using  Kelvin’s  circulation  theorem,  which  states  that  the  circulation  around  a  closed  curve 
moving  with  a  fluid  remains  constant  with  time  for  an  inviscid  barotropic  flow  with  conservative 
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body  forces.  Lifting  bodies  are  represented  as  a  diserete  eollection  of  bound  vortiees  and  an 
equivalent  number  of  control  points.  At  each  moment  in  time,  cireulation  strengths  of  the  bound 
vortiees  are  determined  by  enforcing  flow  tangency  at  the  control  points.  For  a  given  vortex 
fdament,  the  induced  velocity  v  is  given  by  the  Biot-Savart  law  as: 

1  f Tdsxr 

where  F  is  the  circulation  strength  of  the  fdament,  r  is  the  displacement  vector  from  the  filament 
to  the  point  of  interest,  and  ds  is  a  differential  element  of  the  filament.  At  eaeh  iteration  in  time, 
vortiees  are  shed  from  the  lifting  surface  and  eonvected  into  a  wake  so  as  to  preserve  a  net 
cireulation  strength  of  zero  for  the  complete  system. 

We  have  shown  that  the  UVLM  gives  good  agreement  with  data  collected  from  pitehing  and 
plunging  airfoils  in  a  water  tunnel*^.  For  example.  Figure  2  shows  a  eomparison  of  a  plunging 
SD7003  airfoil.  In  this  case,  the  pure  plunging  occurs  at  a  reduced  frequeney  of  k  =  3.93  with  an 
amplitude  of  9.2%  of  the  airfoil  ehord  length.  In  this  view,  the  vertieal  axis  shows  the  solution  at 
a  fixed  location  one  ehord  length  behind  the  trailing  edge  and  the  horizontal  axis  represents  time, 
increasing  to  the  right.  Flow  in  the  tunnel  is  reeorded  using  injected  blue  dye  and  the 
computational  vortices  are  traced  using  colored  dots.  Similar  qualitative  eomparisons  were  found 
for  other  airfoil  motions  and  reduced  frequeneies. 


Figure  2.  Slit  camera  view  comparison  of  UVLM  results  with  water  tunnel  experimental 

data  for  a  pure  plunging  motion. 


3.  Navier-Stokes  Aerodynamics 

Low-order  aerodynamics  models  sueh  as  the  blade  element  and  UVLM  methods  deseribed 
above  make  simplifying  assumptions  that  might  eause  numerical  predictions  to  miss  important 
physical  behavior.  Higher-order  methods  sueh  as  the  Reynolds-averaged  Navier-Stokes  (RANS) 
and  large  eddy  simulation  (LES)  are  appropriate  for  sueh  situations.  Even  where  lower-order 
methods  can  suitably  capture  the  relevant  physics,  Navier-Stokes  methods  can  be  useful  for 
conducting  detailed  inspeetions  of  the  aerodynamic  behavior  in  narrow  regions  of  interest  in  the 
parameter  space. 
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We  have  developed  a  coupled,  aeroelastic  solver  based  on  the  OVERFLOW  Navier-Stokes 
flow  solver  and  a  custom  mode-based  structural  solver'^.  The  solver  uses  an  explicit  coupling 
between  the  flow  solver  and  structural  model.  OVERFLOW  is  a  compressible  Navier-Stokes 
flow  solver  for  structured  overset  solution  domains.  It  is  based  on  a  finite  difference  formulation 
and  is  capable  of  fifth-order  spatial  accuracy  using  the  Weighted  Essentially  Non-Oscillatory 
(WENO)  scheme  and  second-order  temporal  accuracy  using  dual  time  stepping*"^’ For  low 
Mach  number  solutions  or  solutions  encompassing  mixed  Mach  number  flows,  the  code 
incorporates  low  Mach  number  preconditioning'^.  The  implemented  preconditioning  is  time- 
accurate  when  used  with  a  dual  time  stepping  solution.  Multibody  and  moving  body  solution 
capabilities  are  accommodated  using  a  chimera,  or  overset,  solution  approach*^.  Arbitrary  rigid- 
body  motion  may  be  prescribed  and  a  6DOF  model  is  available  for  computing  rigid-body  motion 
due  to  aerodynamic  forces.  The  OVERFLOW  family  of  solvers  has  some  build-in  capability  to 
model  flexible  bodies  and  we  extended  the  program  for  application  to  MAVs  with  flexible  beam¬ 
like  wings. 

The  OVERFLOW-based  aeroelastic  solver  has  been  tested  against  a  two-dimensional  case 
based  on  the  work  of  Heathcote  and  Gursul,  who  examined  the  influence  of  structural  flexibility 
on  thrust  coefficient  and  propulsive  efficiencies  for  a  heaving  airfoil  in  water,  across  a  range  of 
Strouhal  and  Reynolds  numbers.  The  experimental  setup  utilized  an  airfoil  with  an  aluminum 
teardrop  stiffener  along  the  first  third  of  the  chord  and  a  flexible  steel  plate  for  the  last  two-thirds 
of  the  airfoil,  as  shown  in  Figure  3  together  with  a  view  of  near-body  mesh  as  the  body  deforms. 


10% 


I 


H 

3  I  2.5% 


(a)  Airfoil  geometry 


(b)  Deformed  near-body  mesh 


Figure  3.  Navier-Stokes  test  case  configuration  composed  of  a  rigid  airfoil  and  a  flexible 

plate. 


Figure  4  compares  numerical  and  experimental  results  for  the  following  conditions:  Reynolds 
number  Re  =  9000,  Mach  number  Moo  =  0.05,  plate  thickness  b/c  =  0.0056,  and  a  forced 
plunging  of  the  rigid  teardrop  according  to  /i(t)  =  0.05c(l  —  cos((x)t).  Force  histories 

provide  good  agreement  with  experiment.  However,  the  explicit  modal  solutions  overpredict  the 
tip  deflection,  as  represented  by  the  pitching  angle,  defined  as  the  angle  between  the  freestream 
and  a  line  connecting  the  leading  and  trailing  edges  of  the  foil.  This  variation  in  pitching  angle  is 
accompanied  by  an  associated  shift  in  phase  angle  between  the  leading  and  trailing  edge  motion. 
The  linear  mode-based  structural  solver  is  expected  to  overpredict  the  true  nonlinear  response,  so 
using  a  better,  more  accurate  structural  solver  would  likely  eliminate  some  of  the  disagreement 
with  experimental  data. 
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(a)  Pitch  angle 


(b)  Phase  angle 


Figure  4.  Variations  in  predicted  pitch  and  phase  angles  for  changes  in  Strouhal  number 

for  the  Navier-Stokes  test  case. 


Immersed  boundary  methods  have  also  been  used  to  study  problems  that  involve 
flapping The  method  is  well-suited  to  such  flows  because  the  body  does  not  need  to  conform 
to  the  fluid  mesh,  thereby  eliminating  the  need  for  the  remeshing  process  associated  with 
unstructured  solvers  or  the  dynamic  hole  cutting  associated  with  overset  solvers.  In  exchange  for 
a  simplified  handling  of  the  computational  mesh,  special  boundary  conditions  must  be  used  to 
incorporate  the  body  into  the  flow  solution.  We  have  worked  with  Balaras  and  Dong  to 
incorporate  immersed  boundary  codes  into  our  high-order  aerodynamics  tool  suite. 

B.  Nonlinear  Structural  Mechanics 

Small-scale  aircraft,  particularly  flapping-wing  MAVs,  generally  use  very  thin  wings.  These 
might  be  constructed  of  plies  of  composite  material  or  of  an  even  thinner  flexible  membrane 
material.  Beams  or  battens  are  often  incorporated  into  the  wings  to  provide  additional  stiffness. 
As  shown  in  Figure  5,  insect  wings,  with  their  vein  and  cell  structure,  often  inspire  manufactured 
wing  designs.  Structural  deformations  in  such  wings  can  be  large,  especially  for  flapping,  and 
must  be  considered  to  be  nonlinear.  In  some  cases,  the  true  structure  can  be  reduced  to  a  set  of 
equivalent  nonlinear  beams.  Here,  variational  asymptotic  beam  sectional  analysis  can  be  used  to 
facilitate  construction  of  an  equivalent  nonlinear  beam^^’^^.  In  general,  however,  MAV  wings 
should  be  treated  as  nonlinear  plates  or  membranes,  possibly  combined  with  nonlinear  beam 
elements  for  stiffening.  We  have  developed  nonlinear  finite  element  solvers  that  are  based  on  a 
corotational  approximation  of  the  updated  Lagrangian  procedure^"*'^®.  Both  two-dimensional  (2D) 
beam  and  three-dimensional  (3D)  shell  finite  element  formulations  are  available.  In  the 
corotational  approach,  the  total  displacement  of  an  element  is  divided  into  a  component  due  to 
rigid-body  motion  and  a  component  due  to  pure  displacement^’.  The  contribution  of  the  rigid 
motion  is  removed  from  the  element  computations,  leaving  only  the  elastic  deformation.  This 
makes  it  possible  to  consider  problems  with  large  rotations  but  small  strains.  In  three  dimensions, 
a  three-node  triangular  shell  element  is  used  that  combines  an  optimal  membrane  element  (OPT) 
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and  a  discrete  Kirchoff  theory  (DKT)  plate  bending  element^*.  The  effect  of  nonlinear  stress 
stiffening  is  added  to  the  co-rotational  formulation  by  including  a  geometric  stiffness  matrix. 
Limited  verification  and  validation  of  these  solvers  has  been  completed  and  have  shown  good 
agreement  with  theory. 


Source:  Bugboy52.40  (Wikimedia  Commons) 

(a)  Fly  wing 


(b)  Manufactured  wing 


Figure  5.  Comparison  of  an  insect  wing  and  a  manufactured  wing  for  a  flapping  MAV. 


In  Figure  6,  the  3D  nonlinear  beam  model  is  compared  with  results  from  the  commercial 
finite  element  solver  ABAQUS.  A  beam  with  a  90°  elbow  and  a  solid  circular  cross-section  is 
clamped  at  one  end  and  subjected  to  an  impulse  point  load  at  its  comer.  Inertial  and  aerodynamic 
forces  are  zero,  stractural  damping  is  zero,  and  the  force  impulse  lasts  for  the  duration  of  the 
simulation.  The  time  history  of  the  displacement  at  the  free  end  of  the  elbow  compares  well  with 
the  commercial  code  and  demonstrates  the  ability  of  the  solver  to  accurately  capture  large  3D 
bending  and  twisting  dynamics^^. 

The  3D  nonlinear  shell  model  is  compared  with  ABAQUS  in  Figure  7.  For  this  study,  the 
feathering  axis  mns  through  the  mid-chord  of  the  wing  and  the  entire  root  is  clamped,  as 
indicated  on  the  diagram  within  the  figure^"*.  Kinematic  motions  for  this  problem  are  given  by  the 
expressions 


0  =  (l  —  e  ■  (—n/4)sm(o)t) 

ip  =  (^1  —  e-2ooot2^  .  sin((x)t  —  n/2) 

with  a  flapping  frequency  of  u)  =  100  rad/sec  and  no  viscous  damping.  Results  are  given  in 
terms  of  the  normalized  transverse  and  spanwise  displacement  at  point  A,  located  at  the  leading- 
edge  tip  of  the  wing.  Displacements  are  measured  with  respect  to  the  body-attached  coordinate 
system  and,  therefore,  reflect  stmctural  deformations,  but  not  the  prescribed  rigid  body  motions 
due  to  flapping.  The  match  between  the  shell  model  and  ABAQUS  is  very  good.  The  strongly 
nonlinear  transverse  deformations  are  of  the  same  magnitude  as  the  wing  length,  resulting  in  a 
significant  foreshortening  effect  that  is  accurately  captured  by  the  shell  model. 


8 

Approved  for  public  release;  distribution  unlimited 


6 


01““" 


Figure  6.  Structural  beam  model  validation  of  a  clamped  elbow  beam  subjected  to  an 

impulse  load  at  its  corner. 


Figure  7.  Structural  shell  model  validation  of  a  clamped  plate  subjected  to  flapping  motion. 


C.  Kinematics 

Rigid  motion  in  a  flapping- wing  vehicle  can  be  described  in  several  ways.  We  have  relied  on 
a  couple  of  descriptions. 

1.  Berman  and  Wang  Kinematics 

Following  Berman  and  Wang^°,  we  consider  a  flapping  wing  that  is  able  to  rotate  in  three 
degrees  of  freedom  about  a  hinge  at  the  root  of  the  wing.  The  degrees  of  freedom  are  the  azimuth 
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angle  (p,  elevation  angle  6,  and  pitch  angle  77,  as  shown  in  Figure  8.  The  equations  for  these 
angles  are  taken  to  be: 


sin  ^{K^sm(2nft)) 

0  ~  0m  !  \ 

Sin 

6  =  9^  cos(2n ft  +  Og)  +  9q 

tanh(C^  sin(27r/t  +  775)) 


77  =  q. 


tanh(C^) 


+  rio 


(5) 


These  equations  contain  a  total  of  ten  parameters.  The  parameter  may  vary  between  0  for  a 
sinusoidal  waveform  and  1  for  a  triangular  waveform.  The  parameter  varies  between  0  for  a 
sinusoidal  waveform  and  infinity  for  a  step  function. 


Figure  8.  Berman  and  Wang  flapping  wing  kinematics  and  coordinate  systems. 


2.  Spline-Based  Kinematics 

Basis  splines  (B-splines)  can  be  used  to  represent  kinematic  degrees  of  freedom  in  a  very 
general  way.  B-splines  offer  localized  control  of  the  path  and  the  number  of  possible  defining 
parameters  is  essentially  unbounded^*.  Introduce  a  knot  vector  T  with  m  non-decreasing  real 
elements.  Each  kinematic  degree  of  freedom  S(t)  can  be  expressed  as  a  linear  combination  of 
basis  functions: 


m-p-2 

S(t)=  ^  tG  (6) 

i=0 

Here,  p  is  the  order  of  the  spline.  The  Pi  are  called  the  control  points,  and  there  are  tti  —  p  —  1 
such  points.  Their  positions  become  the  parameters  that  describe  the  kinematic  formulation.  The 
B-spline  basis  functions  can  be  defined  using  the  Cox-de  Boor  recursion  formula: 
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£  oil 


<  t  <  ti+1 
otherwise 

T  _|_  -n  ^ 


^i+p 


^i+p+i  f 
i+p+l  “  ^i+l 


M- 


i+l,p- 


-l(t) 


(7) 


D.  Optimization 

Two  types  of  optimization  algorithm  are  used  in  this  effort:  single-objective  optimization  via 
a  gradient-based  method  and  multi-objective  optimization  via  a  genetic  algorithm.  Non-gradient 
methods  require  more  objective  and  constraint  function  evaluations  than  gradient-based  methods. 
Gradient-based  methods  also  have  good  convergence  properties  within  the  radius  of  convergence 
around  an  optimal  point.  Consequently,  gradient-based  methods  are  typically  favored  for 
problems  in  which  the  derivatives  of  the  objective  and  constraints  are  known  analytically  or  can 
be  computed  quickly  and  accurately  by  finite  difference.  Non-gradient-based  methods  can  be 
preferred  for  problems  where  derivatives  are  unavailable  or  difficult  to  compute  and  for 
problems  in  which  the  presence  of  local  minima  confound  a  gradient-based  search  for  a  global 
optimum. 

1.  Gradient-Based  Method 

The  optimization  problem  can  be  stated  formally  as 

m\nii(x) 

X 

s.t. 

g(x)  <  0  (8) 

/i(x)  =  0 
Xf  <  X  <  X-n 

where  gix)  is  the  objective  function,  g{x)  are  the  inequality  constraints,  /i(x)  are  the  equality 
constraints,  and  x  is  the  vector  of  decision  variables. 

The  gradient-based  modified  method  of  feasible  directions  is  used  for  constrained 
optimization  problems,  as  implemented  in  the  CONMIN  and  DOT  optimization  programs^^’^^. 
These  programs  can  also  use  sequential  linear  programming  and  sequential  quadratic 
programming  methods  for  constrained  optimization  and  Broydon-Fletcher-Goldfarb-Shanno 
variable  metric  and  Fletcher-Reeves  conjugate  gradient  methods  for  unconstrained  optimization. 
These  methods  are  iterative  and  require  evaluation  of  the  chosen  objective  function  and  any 
active  constraints  at  each  point  in  the  optimization  search.  Typically,  about  a  dozen  iterations  are 
needed  to  converge  to  an  optimum  solution. 

2.  Genetic  Algorithm 

The  problem  statement  for  a  multi-objective  optimization  is  the  same  as  that  given  in  (8), 
except  that  the  objective  function  ju(x)  becomes  an  n-dimensional  vector  of  functions.  Instead  of 
a  unique  solution,  the  multi-objective  problem  possesses  a  possibly  infinite  set  of  points  that 
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together  define  the  solution.  These  are  called  the  Pareto  points,  and  are  identified  as  points  in  the 
objective  space  jU*  for  which  there  does  not  exist  another  feasible  design  objective  vector  jU  such 
that  jU;  <  jUj*  for  all  i  E  {1,2, ...  n),  and  fXj  <  ju]  for  at  least  one  index  j  E  {1,2, ...  n}. 

The  Non-dominated  Sorting  Genetic  Algorithm-II  (NSGA-II)  is  used  here  to  solve  multi- 
objective  optimization  problems  via  a  non- gradient-based  evolutionary  algorithm.  The  method 
uses  a  fast  non-dominated  sorting  approach  to  search  for  the  Pareto  front^"^.  Elitism  in  the 
selection  process  speeds  the  solution  process  and  also  helps  to  prevent  the  loss  of  good  solutions. 
As  with  many  evolutionary  algorithms,  the  algorithm  identifies  a  large  number  of  alternative 
solutions  lying  on  or  near  the  Pareto-optimal  front  in  a  single  run. 
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APPLICATIONS 


In  this  section,  we  summarize  a  few  applications  of  our  multidisciplinary  optimization 
methodology  to  the  study  of  significant  problems  associated  with  flapping-wing  MAVs.  The 
reader  is  referred  to  the  references  given  in  Appendix  A  for  other  applications  conducted  under 
the  auspices  of  this  laboratory  task. 

A.  Aeroelastic  Shape,  Structure,  and  Kinematics  Design  and  Parameter  Exploration 

The  wings  of  animals  that  fly  are  often  flexible.  As  with  conventional  aircraft,  rigid  wings 
are  simply  too  heavy  for  flight,  and  a  variety  of  numerical  and  experimental  research  has 
indicated  the  benevolent  effects  of  wing  flexibility  for  both  energy  and  lift  considerations, 
outside  of  the  obvious  weight  savings.  Insects  have  no  muscles  embedded  in  their  wings,  and  so 
flexibility  plays  a  passive  role:  adaptive  feathering,  cambering,  and  chordwise  bending,  and 
various  coupling  mechanisms  have  been  found  for  many  species^^.  Birds  and  bats  are  capable  of 
active  morphing,  but  these  creatures  also  use  passive  mechanisms.  Although  the  role  of 
flexibility  in  flying  animals  is  not  fully  understood,  it  almost  certainly  contributes  to  the  extreme 
agility,  maneuverability,  and  stability  of  biological  flyers.  Successful  exploitation  of  such  passive 
flexibility  mechanisms  in  a  mechanical  flapper  requires  an  ability  to  simultaneously  consider 
changes  in  aerodynamic  shape,  flapping  kinematics,  and  structural  layout.  Considering  the 
structure  alone  is  insufficient  due  to  strong  nonlinear  couplings  with  fluid  flow. 

When  designing  a  flapping  aircraft  capable  of  sustained  stationary  hovering,  three  objectives 
must  be  met.  First,  the  wings  must  generate  enough  aerodynamic  lift  to  balance  the  weight  of  the 
vehicle.  Second,  the  mechanical  stresses  that  develop  within  the  wing  structure  must  not  lead  to 
mechanical  failure.  Third,  the  input  power  to  the  system  must  be  minimized.  This  last  point  is 
particularly  important  because  the  energy  budgets  onboard  a  MAV  are  expected  to  be  extremely 
limited.  The  interaction  of  these  objectives  can  be  examined  numerically  by  forming  an 
inexpensive  aeroelastic  optimization  capability  for  power  minimization  under  trim  and  failure 
constraints  ’  .  The  aeroelastic  optimization  used  for  this  effort  combines  Berman  and  Wang 
aerodynamics  and  Berman  and  Wang  kinematics  with  a  nonlinear  beam  model  based  on  a  co- 
rotational  approximation  of  the  updated  Lagrangian  procedure.  These  numerical  models  are 
described  in  earlier  sections. 

The  elastic  wing  is  divided  along  the  span  into  a  number  of  beam  elements  with  rectangular 
cross  sections,  as  shown  in  Figure  9.  The  structural  beam  and  the  aerodynamic  surface  are 
considered  here  to  be  a  unified  object.  In  other  words,  at  a  given  span  station,  the  beam  width 
and  aerodynamic  chord  length  are  identical.  The  chord  length  C;  and  thickness  of  each  element 
is  allowed  to  change  independently.  The  thickness  of  each  element  has  a  strong  effect  on  the 
structural  and  mass  properties  of  the  wing,  but  will  only  slightly  alter  the  aerodynamic  forces 
through  added  mass  terms  in  the  Berman  and  Wang  aerodynamics  model.  Conversely,  the  chord 
length  of  each  element  will  have  a  strong  effect  on  the  aerodynamic  forces,  but  only  have  a 
moderate  impact  on  the  structure.  These  dependencies  are  explicit;  the  coupled  nature  of  the 
aeroelastic  solver  ensures  that  both  chord  and  thickness  distributions  effect  the  aeroelastic 
equations  of  motion  in  an  implicit  manner.  A  direct  method  is  used  to  evaluate  the  sensitivities  of 
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the  objective  and  constraint  equations  to  changes  in  the  design  variables,  and  gradient-based 
optimization  is  used  to  guide  the  optimization  process. 


Figure  9.  Finite  element  discretization  of  a  nonlinear  beam,  with  variable  chord  length  and 

thickness  for  each  element. 


For  a  beam  with  N  finite  elements,  there  are  2N  design  variables  describing  the  shape  and 
structure.  There  are  an  additional  nine  parameters  describing  the  Berman  and  Wang  kinematics, 
so  the  optimization  problem  has  a  total  of  2A1  -I-  9  design  variables: 


X  —  {Ci  Cjy  /ij  •••  h2  (pm  l^cf)  dm  ^0  Vm  Vs  Vo 


(9) 


Baseline,  lower,  and  upper  bounds  for  these  variables  is  given  in  Table  1. 
The  optimization  problem  can  be  formally  stated  as: 


X 

S.t. 

^ave  —  ^ 

KS{tcr,2^  —  0 

«n  =  0 

^min  —  —  ^max 


tl  1,  ... ,  Ngfgpg 

e  =  1, ... ,  Njjy 


(10) 


where  p  is  the  power  required  to  actuate  a  single  wing,  given  by  the  following  equation: 


p  = 


d/fp  dS  p 
f 


dt  dt 


"f  Paero 


(11) 
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Here,  is  the  kinetic  energy  of  the  wing  due  to  both  rigid  body  and  deformational  motions, 
Se  is  the  strain  energy  of  the  beam  due  only  to  deformational  motions,  and  Paero  is 
aerodynamic  power.  It  is  desired  to  minimize  the  peak  power  draw  of  the  system,  which  occurs 
at  a  time  once  the  motion  has  become  periodic.  It  is  further  desired  to  keep  the  elastic 
stresses  that  develop  throughout  the  beam  below  some  critical  failure  stress.  Including  a  stress 
constraint  for  each  beam  element  is  not  desirable.  However,  a  Kreisselmeier-Steinhauser 
function  can  be  computed  over  the  wing  and  a  constraint  applied  to  that  value. 

1 

KS  =  --\n 

K 

Where  is  the  maximum  Von  Mises  stress  within  the  beam  element,  a*  is  the  maximum 
allowable  stress,  and  k  =  20.  The  composite  KS  function  peaks  at  the  critical  time  tcr,2»  which 
time  the  function  must  be  negative.  Unlike  the  power  objective,  though,  tcr,2  can  occur  at  any 
time  during  the  simulation,  not  just  during  the  periodic  flapping  motion.  The  cycle-averaged  lift 
force  L  generated  by  one  wing  must  be  greater  than  some  specified  value  L*.  Also,  all  design 
variables  must  lie  between  side  constraints,  as  given  in  Table  1.  Finally,  the  equations  of  motion 
must  be  satisfied  at  each  time  step  or,  in  other  words,  /?„  =  0  for  each  time  step  n. 


Table  1.  Baseline  values  and  side  constraints  for  the  shape,  structure,  and  kinematics 

design  variables. 


Cl  (mm) 

hi  (mm) 

(pm  (rad) 

9m  (rad) 

9s  (rad) 

minimum 

10 

0.2 

— 7r/2 

0.01 

—njA: 

—njA: 

baseline 

25 

0.6 

nl3 

0.01 

0 

0 

maximum 

60 

2.5 

7r/2 

0.50 

7r/4 

7r/4 

6>o  (rad) 

Vm  (rad) 

rjs  (rad) 

ho  (rad) 

c^(-) 

minimum 

— 7r/2 

— 7r/2 

— 7r/2 

— 7r/2 

0.01 

baseline 

0 

7r/4 

7r/2 

— 7r/2 

0.10 

maximum 

7r/2 

7r/2 

7r/2 

7r/2 

1.00 

The  number  of  beam  elements  is  taken  to  be  Ai  =  10,  which  gives  a  total  of  29  possible 
design  variables.  Each  flapping  cycle  is  broken  into  100  time  steps,  and  8  flapping  cycles  were 
found  to  be  sufficient  to  develop  a  time-periodic  aeroelastic  response.  The  remaining  system 
parameters  are  given  in  Table  2. 
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Table  2.  Parameters  for  the  shape,  structure,  and  kinematics  optimization  of  a  flapping 

wing  in  hover. 


Parameter 

Value 

Parameter 

Value 

N 

10 

Pf 

1.225  kg/m^ 

Wing  length 

10  cm 

Ct 

1.83 

Elastic  modulus 

70  GPa 

Cd(0) 

0.21 

Shear  modulus 

26.9  GPa 

Coin 12) 

3.35 

Wing  density 

2800  kgW 

Cr 

n 

C 

20M 

Ml 

0.2 

a* 

350  MPa 

y-2 

0 

L* 

0.15  N 

f 

19.9  Hz 

Three  optimization  studies  are  considered  and  are  summarized  in  Table  3.  Case  A  considers  only 
the  kinematic  design  variables  and  fixes  the  chord  lengths  and  thicknesses  of  the  beam  elements 
to  a  baseline  design,  leaving  9  active  design  variables.  Case  B  considers  only  the  chord  lengths 
and  thicknesses  that  describe  the  aerodynamic  and  structural  properties  of  the  wing  while 
holding  the  kinematics  to  a  baseline  design.  This  gives  20  active  design  variables.  Finally,  for 
case  C,  all  variables  are  active. 


Table  3.  Cases  for  the  shape,  structure,  and  kinematics  optimization  of  a  flapping  wing  in 

hover. 


Case  N[)y  Description 

A  9  Optimize  kinematics,  fix  aerodynamic  and  structural  geometry 

B  20  Optimize  aerodynamic  and  structural  geometry,  fix  kinematics 

C  29  Optimize  kinematics  as  well  as  aerodynamic  and  structural  geometry 


The  optimal  designs  are  given  in  Figure  10  and  Figure  11.  The  former  gives  the  chord  length 
and  thickness  distribution  along  the  beam  for  cases  B  and  C.  The  shape  for  case  A  is  that  of  the 
baseline.  The  latter  figure  gives  the  kinematic  flapping  motions  for  cases  A  and  C.  The 
kinematics  for  case  B  is  that  of  the  baseline.  Focusing  on  the  shape  and  structure,  both  optimal 
designs  drop  the  thickness  below  the  baseline  value.  The  wing  tip  thickness  for  both  designs  is 
that  of  the  lower  bound  on  thickness,  presumably  to  decrease  the  wing  moment  of  inertia  and 
thus  the  inertial  power.  Both  designs  keep  the  chord  length  of  the  first  few  beam  elements  near 
the  root  close  to  the  baseline  value.  Larger  chord  lengths  are  used  outboard  for  lift  generation, 
where  the  translational  motion  is  larger.  Significant  differences  are  seen  between  cases  B  and  C 
at  the  wing  tip,  where  the  former  pushes  the  length  to  the  lower  bound  on  chord,  while  the  latter 
keeps  the  length  at  a  moderate  value. 
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baseline 


case  B 


case  C 


Figure  10.  Optimal  chord  and  thickness  distributions  for  a  flapping  wing  in  hover. 


baseline  - case  A  - case  C 


t/T  t/T  t/T 

Figure  11.  Optimal  kinematic  flapping  motions  for  a  flapping  wing  in  hover. 


As  opposed  to  the  shape  and  structural  design  variables,  the  optimizer  does  not  significantly 
alter  the  flapping  motions  for  power  reduction.  The  overall  out-of-phase  relationship  between  the 
stroke  motion  0  and  the  wing  feathering  rj  is  preserved,  presumably  for  lift  generation 
requirements.  The  optimal  designs  use  the  stroke  plane  deviation  angle  6  to  push  the  rigid 
flapping  motion  slightly  above  or  below  the  mean  plane,  and  passive  aeroelastic  deformation 
provides  further  deviation.  Case  C  is  able  to  drop  the  amplitude  of  stroke  motion  significantly 
below  the  baseline  value. 

The  time-periodic  power  history  that  results  from  the  optimal  designs  is  shown  in  Figure  12, 
together  with  that  for  the  baseline  design.  Case  A  is  unable  to  decrease  the  peak  power  from  the 
baseline  value  in  a  feasible  manner,  as  the  lift  constraint  must  be  met.  The  kinematic  motions 
alone  were  not  able  to  provide  efficient  solutions.  For  case  B,  fixing  the  kinematic  motions  and 
permitting  changes  to  the  shape  and  structure  provides  a  significant  drop  in  peak  power.  For  case 
C,  also  including  the  kinematic  variables  provides  a  still  greater  reduction  in  peak  power.  A  plot 
of  the  terms  that  make  up  the  total  power  for  case  C  is  instructive.  In  Figure  13,  no  single  term 
dominates,  indicating  that  aerodynamic  and  inertial  forces  are  of  similar  magnitude.  Drops  in  the 
kinetic  energy  rate,  as  the  wing  generates  lift  through  the  mid-stroke,  offset  peaks  in  the 
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aerodynamic  power,  as  energy  is  continuously  stored  and  released  throughout  the  vibrating 
beam,  attenuating  the  total  peak  value. 


Figure  12.  Time-periodic  power  history  for  baseline  and  optimal  designs  of  a  flapping  wing 

in  hover. 


t/T 


Figure  13.  Contributions  to  the  power  for  optimal  case  C. 


The  aeroelastic  wing  motion  of  case  C  is  given  in  Figure  14,  demonstrating  the  ability  of  the 
design  to  leverage  aeroelastic  deformation  in  a  power-optimal  manner.  The  commanded 
kinematics,  as  indicated  by  the  rigid  motion,  provide  a  very  small  amplitude  in  both  stroke 
motion  0  and  stroke  plane  deviation  angle  0,  thereby  decreasing  the  inertial  power  requirement 
associated  with  the  rigid  motion,  dKE/dt.  The  resulting  flexible  motion  substantially  increases 
the  amplitude  of  both  angles,  providing  a  wide  figure-8  motion.  As  noted  earlier,  the  feathering 
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motions  r]  decreases  the  minimum  angle  of  attack  of  the  rigid  wing.  Such  a  rotation  increases  the 
stroke  deviation,  heightening  the  figure-8  for  lift  generation.  This  force  development  increases 
the  aerodynamic  power,  but  judicious  phasing  of  the  elastic  energy  storage  within  the  cycle 
attenuates  the  peak. 


Y,(m)  Y,(m) 

Figure  14.  Rigid  and  flexible  motion  of  the  wingtip  for  case  C. 


In  this  application  of  an  optimization  strategy  to  understand  the  physics  of  flapping  wings, 
three  disparate  types  of  design  variables  are  considered:  aerodynamic  shape,  wing  structure,  and 
kinematics.  The  optimization  shows  that  judicious  tailoring  of  the  parameters  creates  beneficial 
interactions  in  the  complex  coupling  of  aerodynamics  and  structure,  making  it  possible  to  reduce 
power  requirements  without  compromising  either  lift  generation  for  trim  or  structural  integrity. 
The  following  conclusions  can  be  drawn: 

1.  Flexible  wing  motion  can  differ  substantially  from  the  commanded  kinematics  enforced  at 
the  wing  root,  providing  wider  stroke  angles,  higher  stroke  velocities,  and  larger  stroke  plane 
deviations.  Optimal  designs  rely  upon  large  passive  wing  deformation  to  satisfy  lift 
requirements. 

2.  Aerodynamic  power  requirements  are  typically  positive,  while  inertial  and  elastic 
contributions  are  liable  to  be  negative.  The  peak  power  can  be  reduced  by  properly  timing  the 
phase  angle  between  the  aerodynamic  power  and  the  strain  and  kinetic  energy  rates  of  the 
vibrating  beam. 

3.  Changes  to  the  shape  and  structure  are  far  more  effective  at  reducing  power  requirements 
than  changes  to  the  flapping  kinematics. 

These  conclusions  indicate  the  importance  of  an  treating  a  flapping  wing  MAV  as  an  aeroelastic 
system.  An  aerodynamics-only  or  a  structures-only  treatment  will  miss  important  physical 
interactions  that  generate  superior  performance. 

B.  Conceptual  Design  of  Compliant  Mechanisms  for  Flapping  Wings 

Flapping- wing  MAVs  are  generally  of  such  small  scale  and  light  weight  that  the  mechanism 
used  to  actuate  flapping  cannot  be  isolated  from  the  aeroelastic  system.  Aerodynamic  and  inertial 
loads  on  the  wings  are  carried  into  the  actuation  mechanism  through  the  wing  attachment,  where 
they  can  have  a  first-order  effect  on  the  behavior  of  the  mechanism.  For  this  reason,  we  believe 
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that  mechanism  design  should  be  incorporated  into  the  aeroelastic  design  process  for  physical 
exploration.  In  a  preliminary  step  towards  such  an  integrated  capability,  we  considered  the 
conceptual  design  of  the  actuator  alone,  while  subjected  to  aerodynamic,  inertial,  and  structural 
loads  from  attached  flapping  wings^’. 

1.  Optimizing  for  Thrust  Alone 

The  problem  is  illustrated  in  Figure  15.  The  compliant  actuation  mechanism  is  constrained  to 
reside  within  a  2D  box,  shown  in  the  figure  by  the  grey  square.  The  mechanism  is  attached  to 
two  wings  at  the  upper  left  and  right  comers  of  the  box  and  the  wings  are  made  to  flap  by 
applying  a  sinusoidal  point  load  vertically  at  the  bottom  center  of  the  box,  as  indicated  by  the 
applied  force  f.  The  forces  /2  and  shown  in  the  figure  are  fixed  at  zero.  This  produces  a 
simple  flapping  motion  within  the  plane  of  the  actuator  box  that,  if  large  enough,  will  produce  a 
propulsive  force  perpendicular  to  the  box  that  is  sufficient  to  overcome  vehicle  drag. 


Figure  15.  Schematic  of  an  actuated  flapping  wing  in  forward  flight. 


The  actuator  design  can  be  treated  as  a  topological  optimization  problem  in  which  the 
decision  variables  include  both  the  mechanism  topology  and  its  supports.  The  optimization  is 
expressed  here  in  a  gradient-based  form  and  is  solved  using  the  method  of  moving  asymptotes 
(MMA)^^.  The  actuator  region  is  divided  into  a  number  of  elements  and  the  MMA  is  allowed  to 
choose  the  elements  that  will  be  occupied  with  material.  The  problem  can  be  stated  as  follows: 
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min{-Cr,ave) 

X,Z  ^  ^ 

S.t. 

V{x)  <  V* 
5(z)  <  5* 


■  (^1  —  x)  <  grey* 
|u“|  <  u* 

n  =  1,  .. 

■ '  ^steps 

^min  —  1 

e  =  l,. 

..,N, 

^min  —  ■^e  —  1 

e  =  1, .. 

..Ne 

(13) 


The  objective  is  to  maximize  the  cycle-averaged  thrust  coefficient  Cj-ave  while  allowing  the 
elemental  densities  x  and  support  variables  z  to  vary.  The  total  volume  V  of  the  mechanism  is 
required  to  be  below  some  specified  value  and  the  area  of  support  S  is  similarly  required  to  fall 
below  some  value.  It  is  further  required  that  the  magnitude  of  the  displacement  at  the  input 
point  load  be  less  than  some  specified  tolerance  at  each  time  step  to  stabilize  the  optimization 
and  to  provide  a  suitable  load  path  between  the  input  and  the  wing  motion.  An  additional 
constraint  on  the  elemental  densities  essentially  drives  the  elements  to  binary  empty  or  filled 
states.  Lastly,  side  constraints  are  imposed  on  the  decision  variables. 

Symmetry  about  the  x-z  plane  is  assumed,  so  only  half  of  the  model  need  be  modeled.  The 
mechanism  and  the  wing  rod  are  idealized  as  a  single  finite  element  model,  discretized  into 
square  Q4  elements  that  include  geometrically  nonlinear  behavior  to  support  large  elastic 
deformations  of  slender  structures'^.  Aerodynamic  forces  are  applied  to  the  nodes  running  along 
the  lower  surface  of  the  wing  rod.  The  wing  is  assumed  to  be  perfectly  rigid  in  both  a  chordwise 
and  torsional  sense  so  as  to  obtain  a  two-dimensional  response;  only  external  aerodynamic  forces 
in  the  y-z  plane  are  applied  to  the  structure.  The  material  properties  used  for  the  structural  model 
are  shown  in  Table  4.  The  method  of  Buhl  is  used  to  support  the  design  of  the  underlying  support 
structure  along  with  the  topological  design  of  the  mechanism  itself*®.  Within  each  bilinear 
element  that  falls  within  the  actuator  box,  two  springs  are  attached  to  each  node  and  fixed  points 
in  space:  one  spring  aligned  horizontally  and  the  other  vertically.  The  stiffnesses  of  the  eight 
springs  are  taken  to  be  identical,  and  this  value  is  used  as  a  second  design  variable  for  the 
element,  in  addition  to  the  element  density.  The  blade  element  model  of  Berman  and  Wang  is 
used  to  model  the  aerodynamic  forces.  Parameters  describing  the  forward  flight  and  blade 
element  aerodynamic  model  are  shown  in  Table  5.  These  parameters  correspond  to  a  Reynolds 
number  of  Re  =  27,750  and  a  reduced  frequency  based  on  the  half-chord  of  k  =  0.205.  The 
final  equations  of  motion  resulting  from  the  coupling  of  the  structural  and  aerodynamic  models 
form  a  second-order  nonlinear  differential  equation.  This  is  solved  using  an  implicit  time 
marching  scheme  with  Newton-Raphson  subiterations  within  each  time  step  to  drive  the 
nonlinear  residual  below  a  specified  tolerance. 

The  optimization  algorithm  requires  sensitivities  of  the  objective  and  constraint  equations 
with  respect  to  changes  in  the  decision  variables.  A  direct  method  is  used  here  to  evaluate  these 
gradients  analytically.  This  requires  that  a  differential  equation  with  multiple  right-hand  sides 
(one  per  decision  variable)  be  integrated  in  time  along  with  the  system  response.  The  Jacobian 
that  is  needed  by  this  equation  will  be  known  from  the  system  response  computation  at  each  time 
step.  The  Jacobian  is  fairly  sparse  and  well  banded,  so  the  solution  is  not  cumbersome. 
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Optimal  topologies  are  shown  in  Figure  16.  For  eases  1  to  3,  the  underlying  supports,  shown 
in  the  figure  with  blue  eross-hatching,  were  assumed  and  fixed.  In  eases  4  and  5,  the  support 
structure  was  designed  simultaneously  with  the  topology  of  the  actuator.  Within  each  of  these 
problem  types,  the  different  optimal  solutions  represent  different  allowable  volume  fractions.  All 
constraints  were  active  for  each  design.  Table  6  lists  the  key  constraints  and  metrics  for  the 
optimal  designs.  The  constraints  are  the  volume  fraction  constraint  V*  and  the  displacement 
constraint  u*,  with  maximization  of  the  average  thrust  coefficient  Cj  ^ve  the  objective.  From 
this  table,  it  can  be  seen  that  cases  1  and  4  share  the  same  volume  and  displacement  constraints 
and  differ  only  in  the  treatment  of  the  supports.  Cases  3  and  5  are  similarly  paired. 


Table  4.  Structural  parameters  for  the  compliant  actuation  mechanism  design. 


Mechanism 

Wing  rod 

Material 

ABS  plastic 

Carbon  fiber 

Elastic  modulus 

2GPa 

40  GPa 

Poisson’s  ratio 

0.3 

0.3 

Density 

1000  kg/m^ 

1400  kgW 

Thickness 

2  mm 

1  mm 

Failure  stress 

35  MPa  (yield) 

80  MPa  (matrix  cracking) 

Damping, 

5 

5 

Table  5.  Aerodynamic  parameters  for  the  compliant  actuation  mechanism  design. 


Parameter 

Value 

Parameter 

Value 

Chord,  c 

4.1  cm 

Angle  of  attack,  a 

0° 

Wing  length,  1 

12.5  cm 

Actuation  force,  / 

15  N 

Flow  velocity,  I/oo 

10  m/s 

Frequency,  o) 

100  rad/s 

Pf 

1.225  kg/m^ 

Ml 

0.2 

P'XZ 

n 

y-2 

0.2 

Cd(0) 

0.04 

0  kg/m 

Coinn) 

2 

Ay/ 

1.671  ■  10-3 

Cr 

n 

la 

9.062  ■  10-3 
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(a)  Case  1 


(b)  Case  2  (c)  Case  3 


(d)  Case  4  (e)  Case  5 


Figure  16.  Optimal  topologies  for  a  flapping-wing  compliant  actuation  mechanism. 
Hatched  areas  are  clamped.  Cases  1-3  have  fixed  supports;  cases  4  and  5  have  variable 

supports. 


In  case  1,  an  upwards  applied  force  at  the  input  port  creates  a  positive  moment  about  the  wing 
joint,  leading  to  a  negative  flapping  angle  and  a  downward  motion  of  the  wing  tip.  This  behavior 
is  also  observed  in  the  corresponding  variable  support  design  case  4.  Cases  2,  3,  and  5,  which 
have  larger  volume  fractions,  exhibit  opposite  behavior:  an  upward  force  produces  a  positive 
flapping  angle. 

The  variable  support  designs  are  simpler  than  the  fixed  support  designs,  with  fewer  moving 
parts  and  superior  performance  to  the  corresponding  fixed  support  design.  The  variable  support 
design  of  case  4  is,  for  example,  exhibits  roughly  four  times  more  thrust  than  the  corresponding 
fixed  support  design  of  case  1.  It  is  not  clear  what  purpose  the  truss-like  structure  of  case  5 
serves,  as  it  does  not  bear  any  stresses  and  does  not  appear  to  contribute  to  the  kinematics. 
Potentially,  this  structure  is  an  artifact  of  the  optimization  process  and  appears  to  fill  out  the 
volume  constraint  without  adversely  impacting  the  compliance  of  the  mechanism. 


Table  6.  Summary  of  performance  metrics  for  optimal  topologies  for  a  flapping- wing 

compliant  actuation  mechanism 


Constraints 

Performance 

Case 

7* 

u*  (mm) 

^L,ave 

^T,ave 

^P,ave 

^mech  (MPa) 

^wing  (MPa) 

1 

0.170 

2 

0.0019 

0.0178 

0.0687 

25.85 

18.79 

2 

0.235 

2 

0.0042 

0.0977 

0.1129 

8.68 

27.62 

3 

0.300 

2 

0.0030 

0.0735 

0.0989 

10.45 

24.80 

4 

0.170 

2 

0.0015 

0.0716 

0.0980 

11.65 

26.58 

5 

0.300 

2 

0.0028 

0.0909 

0.1098 

14.27 

32.78 

Changes  in  the  behavior  of  the  system  for  off-design  actuation  are  shown  in  Figure  17  and 
Figure  18.  For  a  fixed  actuation  frequency  of  100  rad/sec.  The  left-hand  plot  in  Figure  17  shows 
the  effect  of  changes  in  the  input  force  on  the  flapping  amplitude.  The  design  point  of  15  N  is 
shown  in  the  figure  by  a  dashed  line.  A  monolithic  relationship  exists  between  the  force  and 
flapping  rotation,  with  geometric  nonlinearities  and/or  aerodynamic  damping  hardening  the 
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response  under  higher  input  forces.  The  right-hand  plot  in  the  figure  shows  the  peak  von  Mises 
stress  within  the  mechanism  as  a  function  of  the  cycle-averaged  thrust  coefficient.  Each  trend  of 
this  plot  starts  with  zero  input  actuation,  and  so  the  only  force  present  along  the  fixed  flat-plate 
wing  at  zero  angle  of  attack  is  the  viscous  drag.  As  such,  each  case  starts  with  a  thrust  coefficient 
of —0.04,  which  is  the  specified  CoCO)  from  Table  5.  The  flapping  motion  initiates  a  thrust  force 
that  counteracts  the  drag.  Each  case  is  able  to  produce  a  net  positive  thrust  except  case  1 .  From 
the  vantage  point  of  Figure  17,  case  2  is  superior  to  the  remaining  designs  in  that  it  provides  the 
highest  level  of  propulsion  with  the  lowest  mechanism  stresses,  which  should  translate  to  a  long 
fatigue  life. 

Figure  18  shows  results  of  fixing  the  input  actuation  at  15  N  and  traces  the  response  as  a 
function  of  the  actuation  frequency.  As  this  frequency  approaches  zero,  the  motion  is  quasi¬ 
steady  and  only  the  input  and  elastic  forces  contribute.  The  resulting  deformation  is  large  and 
nonlinear.  In  particular,  case  2  flaps  at  an  amplitude  of  75°.  For  nonzero  frequencies,  the 
aerodynamic  forces  resist  the  motion  and  the  inertial  wing  loads  require  a  greater  power  input. 
The  net  result  is  a  decreased  effectiveness  of  the  compliant  mechanism  with  the  introduction  of 
aerodynamic  and  inertial  loads,  which  has  also  been  noted  experimentally.  The  thrust  data  on  the 
right  of  the  figure  begins  at  (0)  for  the  same  reasons  as  above,  and  sees  an  initial  quadratic 
growth  with  frequency.  This  growth  is  diminished,  and  eventually  reversed,  by  the  loss  of 
flapping  amplitude.  It  should  be  noted  that  this  is  strictly  a  coupled  aeroelastic  phenomena, 
which  has  also  been  observed  experimentally;  increasing  the  actuation  frequency  in  a  vacuum 
would  result  in  a  quadratic  growth  in  thrust,  eventually  curbed  by  excessive  geometric 
nonlinearities  or  inertial  loads.  From  the  vantage  point  of  Figure  18,  case  5  is  superior  and  is  able 
to  achieve  a  high  peak  thrust  in  the  high-frequency  domain,  although  case  2  is  the  best  at  lower 
frequencies.  The  stresses  for  case  6  (not  shown)  are  relatively  low.  Both  effects  are  thought  to  be 
due  to  the  rigid  body-like  motions  emulated  by  this  variable  support  topology,  which  is  less 
susceptible  to  aerodynamic  resistance  and  elastic  failure. 

The  inclusion  of  geometrically  nonlinear  elastic  forces,  unsteady  inertial  forces,  quasi-steady 
aerodynamic  forces,  and  oscillating  actuation  forces  into  a  gradient-based  topology  optimization 
scheme  has  produced  compliant  actuation  mechanism  designs  that  provide  required  levels  of 
thrust  while  including  the  coupling  that  exists  between  the  flapping  wing  and  the  mechanism 
itself  Modifying  the  constraint  boundaries  and  the  support  details  provides  greater  control  over 
the  stresses  than  topological  optimization  alone.  This  is  an  important  consideration  from  a 
fatigue  standpoint,  where  cyclic  high-speed  stresses  can  greatly  reduce  the  expected  life  and 
reliability  of  a  flapping  aircraft.  Variable  supports,  in  particular,  were  shown  to  be  effective  for 
creating  simple  designs  with  few  moving  parts. 
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■e —  csist;  I  — B —  case  2  — ^ —  ease  1  O  ease  4  — • —  case  5 


Figure  17.  Variation  in  flapping  performance  due  to  changes  in  compliant  actuator  input 

force  (cD  =  100  rad/sec). 


■© —  asst;  I  — B —  ease  2  — * —  ease  1  O  ease  4  — • —  case  5 


Figure  18.  Variation  in  flapping  performance  due  to  changes  in  compliant  actuator  input 

frequency  (/"=  15  N). 


2.  Optimizing  for  Thrust  and  Rolling  Moment 

The  work  described  above  assumes  symmetric  flapping  motion  and  aims  to  maximize  thrust 
production  for  efficient,  sustained,  steady,  forward  flight.  This  work  has  been  extended  to 
consider  the  optimization  of  a  compliant  mechanism  designed  to  perform  two  tasks"^^.  First,  as 
before,  the  mechanism  must  be  able  to  transform  an  input  actuation  force  into  symmetric 
flapping  motions  of  left  and  right  wings.  The  flapping  amplitude  must  be  large  enough  to 
produce  sufficient  thrust  to  overcome  the  aerodynamic  drag  and  achieve  sustained  forward  flight. 
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Second,  under  different  actuation,  the  mechanism  must  be  able  to  produce  unsymmetric  motion 
between  the  left  and  right  wings  that  gives  rise  to  an  aerodynamic  rolling  moment.  This  will 
allow  turning  maneuvers.  The  topology  of  the  mechanism  itself  must  remain  symmetric  to 
facilitate  the  first  goal,  as  well  as  to  allow  for  equal  capacity  in  left  and  right  turns.  These  two 
goals  are  viewed  as  multiple  load  cases  and  are  handled  as  a  multi-objective  optimization 
problem. 

The  geometry  of  the  problem  is  that  of  Figure  15,  but  now  the  entire  configuration  must  be 
modeled  so  as  to  include  unsymmetric  actuation  and  response.  Input  actuation  forces  may  be 
applied  according  to  one  of  two  scenarios.  In  the  first  scenario,  separate  actuators  are  utilized  for 
the  symmetric  and  unsymmetric  flapping  motion.  A  single  oscillating  point  load  acting  along 
the  centerline  is  used  to  generate  symmetric  motion.  A  single  oscillating  point  load  f^,  offset 
from  the  centerline,  is  used  to  produce  unsymmetric  motion.  In  the  second  scenario,  equal 
oscillating  point  loads  at  /2  and  are  used  to  produce  symmetric  flapping  and  a  point  load  at 
is  again  used  to  produce  unsymmetric  flapping.  The  two  scenarios  are  shown  graphically  in 
Figure  19. 


Figure  19.  Two  scenarios  for  symmetric  and  unsymmetric  actuation  of  left  and  right  wings 

using  a  compliant  mechanism. 


The  optimization  problem  is  slightly  modified  from  the  earlier  statement: 

X  ^  ^ 

S.t. 

r  >  r* 

^T,ave  ^ 

V(_x)  <  V*  (14) 

■  (1  —  x)  <  grey* 

I  <  It  TL  1,  ... ,  Aig^gpg 

^min  —  —  f  ^  1/  ■■■  >  Ng 
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Here  it  is  desired  to  maximize  the  peak  rolling  moment  coefficient  Q  during  the  flapping  cycle 
such  that  the  cycle-averaged  thrust  coefficient  is  greater  than  or  equal  to  some  trim  requirement. 
Taking  Cj  =  0  indicates  that  the  propulsive  forces  exactly  cancel  with  the  viscous  drag.  Also,  we 
do  not  here  consider  the  design  of  the  support  locations.  The  structural  and  aerodynamic 
parameters  are  those  from  Table  4  and  Table  5,  except  that  Co(0)  =  0.025. 

Two  cases  are  considered  for  each  actuation  scenario:  one  mechanism  designed  to  a 
maximum  volume  fraction  of  V*  =  0.21  and  the  other  designed  to  a  larger  volume  fraction  of 
V*  =  0.28.  The  maximum  rolling  moment  coefficients  and  cycle-averaged  thrust  coefficients 
obtained  from  the  four  optimal  mechanisms  are  reported  in  Table  7.  The  mechanisms  themselves 
are  shown  in  Figure  20.  In  all  designs,  a  positive  input  actuation  force  targeting  symmetric 
motion  produces  an  upward  flapping  motion  in  the  wings,  as  with  cases  2,  3,  and  5  in  the  thrust- 
only  optimization.  For  unsymmetric  motion,  a  positive  force  on  the  right  side  of  the  mechanism 
produces  an  upward  motion  in  the  right  wing.  How  the  left  wing  reacts  to  this  force  governs  the 
rolling  moment  performance  of  the  mechanism,  with  larger  unsymmetric  behavior  leading  to 
larger  rolling  moments.  In  mechanism  1,  the  left  wing  does  not  respond  significantly  to  actuation 
on  the  right  side:  the  cross  bracing  between  the  two  sides  of  the  mechanism  is  too  stiff  to 
transmit  the  force  across  the  mechanism.  This  can  also  be  seen  in  Figure  21,  which  shows  the 
relationship  between  the  sinusoidal  input  loading  and  the  flapping  angle.  For  mechanism  1,  there 
is  a  90°  phase  lag  between  the  input  force  and  the  wing  motion.  Mechanism  2  is  observed  to  be 
stiff  due  to  the  higher  volume  fraction  constraint,  which  drops  the  amplitude  of  the  flapping 
motions  and  the  associated  aerodynamic  flight  metrics.  The  rolling  moment,  in  particular,  is  less 
than  that  from  mechanism  1,  despite  the  fact  that  the  wing  motions  are  more  unsymmetric.  Also, 
it  can  be  seen  that  the  phase  lag  between  the  input  and  output  is  lower  than  that  in  mechanism  1, 
as  the  aerodynamic  resistance  is  weaker.  Mechanism  3,  in  which  the  symmetric  and  unsymmetric 
actuators  are  shared,  provides  the  largest  rolling  moment  coefficient  of  any  design.  However,  the 
design  is  not  successful  at  thrust  production,  with  a  lower  cycle-averaged  thrust  coefficient  than 
every  other  mechanism.  This  may  be  due  to  the  fact  that  the  topology  has  less  connectivity 
between  the  input  ports  and  the  clamped  boundary  conditions  than  the  other  mechanisms  and, 
thus,  struggles  to  use  these  fixed  points  to  convert  the  input  forces  into  bending  moments  about 
the  wing  joints.  This  indicates  a  topological  tradeoff  between  effectiveness  in  symmetric  and 
unsymmetric  flight.  Mechanism  4,  with  a  higher  volume  fraction  constraint,  has  a  better 
symmetric  thrust  coefficient  than  mechanism  3,  but  is  too  stiff  in  the  unsymmetric  sense,  leading 
to  a  very  small  maximum  rolling  moment. 

The  mechanisms  designed  under  scenario  2  are  topologically  simpler,  but  are  less  successful 
in  generating  aerodynamic  thrust  and  rolling  moments.  From  the  multi-objective  standpoint  of 
being  able  to  perform  both  tasks  moderately  well,  mechanism  1  may  be  considered  the  superior 
design.  From  a  manufacturing  vantage,  however,  mechanism  3  may  be  preferred. 

Several  conclusions  can  be  drawn  from  the  results  of  the  multi-objective  compliant 
mechanism  design  study.  First,  mechanisms  designed  under  high  volume  fraction  constraints  are 
too  stiff  and  are  unable  to  leverage  deformations  to  produce  substantial  aerodynamic  loads. 
Second,  symmetric  actuation  is  able  to  produce  high  thrust  via  large  flapping  amplitudes. 
Unsymmetric  actuation,  on  the  other  hand,  is  more  dependent  on  large  discrepancies  between  the 
right  and  left  wing  motions  to  produce  high  rolling  moments.  Third,  using  the  same  actuators  for 
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both  symmetric  and  unsymmetric  motion  leads  to  topologically  simpler  designs,  but  are  not  as 
successful  in  robustly  producing  the  desired  aerodynamic  loads  as  designs  using  separate 
actuators.  Finally,  successful  operation  under  both  symmetric  and  unsymmetric  load  cases 
suggests  that  it  may  be  possible  to  smoothly  control  shifts  between  longitudinal  and  lateral  flight. 


Table  7.  Summary  of  optimal  compliant  actuation  mechanisms  designed  for  maximum 
rolling  moment  and  a  trimming  cycle-averaged  thrust  coefficient. 


Mechanism 

Scenario 

1/* 

Symmetric  flapping 

r  r 

^L,max  ^T,ave 

Unsymmetric  flapping 

r  r 

^L,max  ^T,ave 

1 

1 

0.21 

~0 

0.0377 

0.4891 

0.0155 

2 

1 

0.28 

~0 

0.0047 

0.3577 

-0.0021 

3 

2 

0.21 

~0 

-0.0033 

0.5224 

0.0312 

4 

2 

0.28 

~0 

0.0028 

0.0642 

0.0034 

mechanism  1 

mechanism  2 

mechanism  3 

mechanism  4 

Figure  20.  Optimal  compliant  actuation  mechanisms  designed  for  maximum  rolling 
moment  and  a  trimming  cycle-averaged  thrust  coefficient. 


mechanism  1 


mechanism  2 


mechanism  3 


mechanism  4 


-10  0  10 
f-sin(o)-t)  (N) 


■  symmetric 


un-symmetric:  right  wing 


un-symmetric:  left  wing 


Figure  21.  Phase  diagrams  of  flapping  rotation  due  to  input  actuation  (positive  angles 

correspond  to  wing  tips  up). 
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C.  Influence  of  Wing  Stiffener  Arrangements  on  Flapping  Wing  Performance 

Inspired  by  nature,  engineers  often  design  mechanical  flappers  with  wings  composed  of  a 
thin,  flexible  membrane  affixed  to  an  underlying  lattice  of  stiffeners.  This  type  of  wing  is 
desirable  for  meeting  stringent  payload  and  power  budget  requirements.  Careful  selection  of  the 
stiffener  pattern  can  provide  additional  aeroelastic  benefits.  Previous  work  in  this  area  includes  a 
description  of  the  aerodynamic  benefit  of  one  skeletal  arrangement  over  another'*'', 
characterization  of  the  benefits  of  perimeter  reinforcement  for  hovering  wings'*^,  and  studies  into 
the  effects  of  variable  spar  stiffness  distribution  in  a  given  wing  topology'**’''*^.  These  research 
efforts  used  experimental  testing,  but  a  comprehensive  experimental  study  detailing  the 
relationship  between  stiffener  arrangements  and  flapping  flight  performance  is  probably 
impractical  given  the  quantity  and  diversity  of  topologies  that  should  be  considered.  A  numerical 
optimization  scheme  offers  a  means  to  locate  preferred  membrane  wing  stiffener  arrangements 
for  a  series  of  aerodynamic  metrics,  as  well  as  identify  the  tradeoff  curves  that  connect  the 
disparate  design  objectives. 

The  conceptual  design  of  the  stiffener  arrangement  in  a  membrane  wing  is  essentially  a 
topological  optimization  problem.  However,  a  solid  isotropic  material  with  penalization  (SIMP) 
method  such  as  that  used  above  for  compliant  mechanism  design  is  considered  undesirable  for 
this  application  because  of  the  potential  existence  of  multiple  local  minima  and  the  need  for 
difficult  and  costly  gradient  computations.  Instead,  the  NSGA-II  genetic  algorithm  is  used 
because  it  offers  gradient-free  global  multi-objective  optimization.  Three  objectives  are 
considered:  maximizing  cycle-averaged  lift,  maximizing  cycle-averaged  thrust,  and  minimizing 
the  cycle-averaged  power  requirement. 

The  genetic  algorithm  acts  on  a  genome  that  encodes  a  developmental  program.  When  that 
program  is  compiled  and  executed,  it  uses  a  Lindenmayer  system  (L-system,  for  short)  to 
develop  a  wing  topology  consisting  of  a  set  of  connected  closed  membrane  cells'*^.  The  edges  of 
the  cells  are  taken  to  be  the  structural  stiffeners.  By  encoding  a  developmental  program  into  the 
genome  instead  of  a  fine  grid  of  potential  stiffeners,  the  number  of  design  variables  may  be  kept 
relatively  small  while  the  attainable  topologies  can  be  extremely  diverse.  Also,  each  design 
variable  has  a  global  influence  on  the  resulting  topology  and  aeroelastic  response,  which  makes 
crossover  and  mutation  operations  in  the  genetic  algorithm  more  efficient. 

An  L-system  is  a  parallel  rewriting  system,  often  used  to  model  growth  processes  in  plants^**. 
Here,  Binary  Propagating  Map  OL-systems  with  markers  (mBPMOL-systems,  for  short)  are  used 
to  define  the  membrane  wing  stiffener  arrangement  and  to  interface  with  the  genetic  algorithm 
for  topology  optimization.  These  are  a  subclass  of  L-systems.  In  a  map  L-system,  the  starting 
seed  must  form  a  closed  cell.  The  method  is  further  termed  “binary”  because  the  cells  are 
required  to  divide  in  two  during  the  division  process,  and  “propagating”  because  cells  cannot 
fuse  together  or  disappear.  The  system  is  defined  by  a  tuple  G  =  (V,  o),  P),  composed  of  an 
alphabet  V  that  consists  of  a  set  of  symbols,  an  axiom  o)  that  defines  the  initial  state  of  the 
system,  and  a  set  of  production  rules  P  that  define  how  variables  change  over  time.  Each 
production  rule  consists  of  predecessor  and  a  successor  strings,  and  these  rules  are  applied 
simultaneously,  once  per  iteration. 

An  example  set  of  production  rules  is: 
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(15) 


A  ^  B[+A]x[-A]B 
B  ^A 

X  ^  X 

Where  the  alphabet  V  =  {A  B  A  production  rule  is  assigned  to  each  symbol  in  the 
alphabet,  with  an  identity  rule  for  x.  Each  letter  outside  the  brackets  specifies  an  edge 
subdivision,  and  the  letter  inside  the  brackets  specifies  a  marker  label  for  a  possible  cell  division. 
The  marker  is  placed  to  the  left  of  the  predecessor  if  a  positive  sign  precedes  the  letter,  and  to  the 
right  if  a  negative  sign  precedes  it.  After  markers  have  been  placed,  each  cell  is  scanned  for 
matching  markers.  When  matching  markers  are  found  on  different  edges  in  a  cell,  they  are  joined 
to  subdivide  the  cell  and  form  a  new  edge.  If  more  than  one  pair  is  found,  the  first  pair  found  that 
creates  angles  and  an  area  larger  than  some  prescribed  lower  bound  is  used.  Remaining  markers 
are  discarded.  Starting  from  an  axiom  o)  =  ABAB,  Figure  22  shows  the  first  four  iterations  of  the 
method. 


initial  map 


first  step 


second  step 


Figure  22.  First  four  iterations  of  an  example  cellular  division  process. 


The  process  can  be  terminated  after  each  cell  is  smaller  than  some  required  area,  after  a  defined 
maximum  number  of  steps,  or  when  the  system  becomes  invariant  with  successive  iterations. 

The  above  cellular  division  method  is  used  in  conjunction  with  the  genetic  algorithm  by 
encoding  it  into  a  vector  of  real  numbers,  which  is  used  as  the  genome.  The  genome  is 
partitioned  as: 


Y  =  (f" 


yo) 


(16) 


30 

Approved  for  public  release;  distribution  unlimited 


where  {Y-^  Y2  •••  encodes  the  axiom,  [y[^  Y^^  •••  encodes  the  first 

production  rule,  and  so  on.  Kobayashi  and  others  provide  methods  by  which  symbols  and 
production  rules  can  be  converted  into  real  numbers^*. 

Once  the  final  cellular  division  topology  has  been  obtained,  the  square  geometry  is 
transformed  into  the  space  defined  by  the  wing  planform.  The  four  outer  cell  walls  become  the 
leading  edge,  wing  root,  trailing  edge,  and  wing  tip.  The  edge  vertices  within  the  four  outer  walls 
are  mapped  into  the  wing  planform,  and  the  edges  are  redrawn  to  keep  the  stiffeners  straight. 

The  wing  is  modeled  structurally  as  a  collection  of  triangular  membrane  elements  with 
embedded  Euler-Bemoulli  beam  elements.  A  typical  wing  is  shown  in  Figure  23. 


Figure  23.  Sample  wing  finite  element  model  for  studying  the  effect  of  stiffener 
arrangements  on  aeroelastic  performance.  Membrane  elements  are  outlined  with  thin 
black  lines  and  beam  elements  are  shown  as  thick  blue  and  red  lines. 


This  model  can  be  considered  accurate  if  the  stresses  that  develop  within  the  membrane  as  a 
result  of  the  loading  are  smaller  than  the  pre-stresses.  This  condition  has  been  demonstrated  by 
experimental  validation  for  membrane  wing  skeletons^^.  Excessive  deformation  may  require 
nonlinear  membrane  models,  but  this  is  not  considered  here.  To  reduce  the  cost  of  the  system 
analysis,  the  solution  vector  is  approximated  as  a  linear  combination  of  a  relatively  small  number 
of  natural  vibration  modes.  The  resulting  structural  equations  of  motion  are  of  the  form 


Mr  ■  rj  +  Cj.  ■  TJ  +  Kj.  ■  TJ  = 


(17) 


The  airloads  are  computed  using  Peters  aerodynamics,  which  allows  for  coupling  with  a 
generic  unsteady  wake  model  and  relies  upon  Glauert  expansions  to  form  matrix  expressions 
analogous  to  Equation  17^^.  The  method  allows  for  large  rigid  body  motions  and  small  elastic 
deformations.  The  vehicle  is  inclined  at  fixed  angle  of  attack  relative  to  the  incoming  flow  and 
the  wings  flap  according  to 
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(18) 


oo 


h(t,  ^ 


hn  ■  cos(n(p) 


n=0 


where  (p  =  cos  ^(x/b)  and  x  is  in  a  local,  body-attached  reference  frame.  Peter’s  inflow  theory 
is  used  to  calculate  the  flow  induced  by  the  trailing  wake^"*. 

The  coupled  aerodynamic  and  structural  equations  are  solved  using  a  spectral  element 
method  that  exploits  the  cyclic  nature  of  the  problem.  The  set  of  differential  equations  describing 
the  aeroelastic  system  are  converted  into  an  algebraic  set  of  monolithic  time-periodic  equations. 
This  method  is  efficient  and  bypasses  the  initial  transients  of  the  system^^’^*’.  It  is  also  possible  to 
use  /jp-refmement  for  temporal  accuracy  and  asymptotic  annihilation^^.  The  main  drawback  is 
that  the  final  algebraic  system  of  equations  may  be  very  large,  equal  to  the  product  of  the  number 
of  spatial  and  temporal  degrees  of  freedom.  Use  of  modal  reduction  for  the  structure  helps 
alleviate  this  cost. 

A  Zimmerman  planform  with  a  root  chord  of  16  cm,  a  tip  chord  of  4  cm,  and  a  wing  length 
(root  to  tip)  of  40  cm  is  used  for  the  following  results.  This  planform  is  shown  in  Figure  23.  A 
constant  parabolic  arc  airfoil  section  is  used,  with  a  peak  camber  of  2%  of  the  local  chord,  no 
twist,  and  no  dihedral.  The  wing  is  oriented  at  an  angle  of  attack  of  a  =  4°,  the  flow  velocity  is 
Uco  =  10  m/s,  and  the  amplitude  of  the  sinusoidal  flapping  motion  is  ^  =  30°  at  a  frequency  of 
o)  =  40  rad/s.  The  resulting  reduced  frequency  is  k  =  0.32  and  the  Reynolds  number  is 
Re  =  10,000.  The  drag  coefficients  are  set  to  C^q  =  0.05  ad  =  2.0. 

The  wing  is  constructed  from  a  thin  latex  membrane  sheet  adhered  to  a  carbon  fiber  skeleton. 
The  material  properties  are  summarized  in  Table  8.  The  membrane  is  assumed  to  have  an 
isotropic  spatially  uniform  pre-tension  and  the  skeleton  is  composed  of  a  lattice  of  unidirectional 
laminate  strips  with  rectangular  cross  sections.  In  order  to  investigate  the  effect  of  skeletal 
stiffness  upon  the  optimal  topological  configuration,  two  cases  are  studied.  The  beam  elements 
that  make  up  the  first  quarter-chord  of  the  root  and  the  first  quarter-chord  of  the  leading  edge  are 
clamped  in  the  body-attached  coordinate  system. 

Three  objective  functions  are  considered:  maximizing  the  cycle-averaged  lift,  maximizing 
the  cycle-averaged  thrust,  and  minimizing  the  cycle-averaged  power  consumption.  The  power  is 
a  sum  of  aerodynamic  power,  strain  energy  rate,  and  kinetic  energy  rate,  though  the  first  is  most 
important  for  forward  flapping  flight,  and  the  only  one  of  the  three  that  has  a  non-zero  average 
over  a  cycle.  No  constraints  are  considered  in  the  optimization.  The  procedure  used  to  find  the 
Pareto  front  spanning  these  three  objectives  is  as  follows.  First,  single-objective  optimization 
runs  are  conducted  with  each  objective.  The  three  optimum  designs  are  then  combined  with  a 
population  of  several  random  designs  to  locate  the  three  Pareto  fronts  spanned  by  the  three  dual¬ 
objective  optimization  problems.  Finally,  these  fronts  are  combined  with  additional  random 
designs  to  locate  the  three-objective  Pareto  front.  A  population  size  of  100  is  used  for  the  first 
two  steps,  and  300  is  used  for  the  final  step,  each  with  a  maximum  of  100  generations. 
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Table  8.  Material  and  geometric  properties  of  the  membrane  wing  and  stiffeners. 


Membrane 

Flexible  case 

Stiff  case 

LE 

Battens 

LE 

Battens 

Elastic  modulus,  E  (MPa) 

2 

3x10^ 

3x10^ 

3x10^ 

3x10^ 

Poisson’s  ratio,  v 

0.5 

0.34 

0.34 

0.34 

0.34 

Density,  p  (kg/m  ) 

1200 

1600 

1600 

1600 

1600 

Thickness  (mm) 

0.1 

1.4 

0.4 

2.0 

0.8 

Width  (mm) 

— 

5.0 

1.0 

5.0 

3.0 

Pre-stress,  N^,  Ny  (N/m) 

10 

— 

— 

— 

— 

The  resulting  extremes  of  the  final  Pareto  front  for  the  flexible  skeleton  parameters  are  given 
in  Figure  24,  corresponding  to  the  highest  average  thrust,  the  lowest  average  power,  and  the 
highest  average  lift.  The  relevant  aerodynamic  metrics  are  provided  for  each,  as  well  as  the 
cellular  representation  and  the  resulting  wing  topology.  The  propulsive  efficiency  rj  is  also  given, 
but  is  not  used  during  the  optimization.  The  topology  of  the  peak  thrust  design  is  oriented  for  the 
wing  to  behave  in  a  beam-like  fashion.  In  contrast,  the  minimum  power  design  is  very  sparse  and 
flexible  in  a  chordwise  bending.  The  peak  lift  design  is  a  combination  of  the  two,  with  a  distinct 
topological  change  at  midspan. 


Figure  24.  Single-objective  optimal  designs  for  a  flexible  wing  skeleton. 


The  optimal  topologies  for  the  stiffer  skeletal  parameters  are  given  in  Figure  25.  The  thrust 
and  power  optimal  designs  are  very  similar  to  the  flexible  case.  The  lift-optimal  design  is, 
however,  very  different,  with  a  large  region  of  free  membrane  skin  surrounded  by  battens.  These 
differences  indicate  that  cross-sectional  properties  should  be  included  in  the  optimization  along 
with  the  topology  variables,  so  as  to  frilly  exploit  the  aeroelastic  behavior. 
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Figure  25.  Single-objective  optimal  designs  for  a  stiff  wing  skeleton. 


Figure  26  and  Figure  27  show  the  thrust-power  Pareto  fronts  for  both  the  flexible  and  stiff 
cases,  respectively.  Several  topologies  are  given  along  the  fronts,  showing  the  gradual  transition 
from  diagonally  oriented  minimum  power  design  to  spanwise  maximum  thrust  design. 
Uniformly  distributed  random  designs  are  also  plotted.  In  the  flexible  case,  all  of  the  random 
designs  are  densely  clustered  just  behind  the  Pareto  front,  indicating  that  off-optimal  designs 
perform  comparably  to  optimal  designs.  Contrastingly,  in  the  stiff  case,  many  of  the  designs  are 
well  removed  from  the  front,  in  general,  and  the  single-objective  optimal  designs,  in  particular. 
This  indicates  that  careful  selection  of  the  wing  stiffener  arrangement  can  have  a  significant 
impact  on  performance  when  the  support  skeleton  contributes  substantially  to  the  total  mass  and 
stiffness  of  the  wing  and,  therefore,  has  a  significant  influence  on  aeroelastic  behavior  of  the 
wing. 


Figure  26.  Thrust-power  Pareto  front  for  the  flexible  case. 
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Figure  27.  Thrust-power  Pareto  front  for  the  stiff  case. 


Since  the  wing  topologies  with  stiffer  battens  have  a  more  diverse  design  space  and  greater 
need  for  optimization,  only  this  configuration  will  be  considered  further.  The  three-objective 
Pareto  front  is  shown  in  Figure  28.  Projections  of  the  front  into  the  lift-thrust  and  thrust-power 
planes  are  also  given.  It  can  be  observed  in  this  figure  that  the  tradeoff  between  lift  and  thrust  is 
not  as  clearly  defined  as  that  between  power  and  thrust.  Designs  with  high  lift  tend  to  also  have 
high  thrust,  and  so  the  conflict  between  these  two  objectives  is  not  strong.  Power  requirements 
for  flapping  can,  however,  typically  be  lowered  by  reducing  the  aerodynamic  resistance  as  much 
as  possible,  which  clearly  conflicts  with  the  generation  of  lift  and  thrust. 


Figure  28.  Thrust-power-lift  Pareto  front  for  the  stiff  configuration. 


The  aeroelastic  response  of  the  thrust-optimal  case  is  shown  in  Figure  29,  which  includes 
nine  snapshots  of  the  wing  throughout  the  flapping  stroke.  One  side  shows  a  translucent  rigid 
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wing  and  a  contoured  flexible  wing,  where  the  contour  is  the  out-of-plane  displacement.  The 
other  side  shows  the  flexible  wing  stiffener  arrangement.  The  deformation  patterns  for  this  case 
are  dominated  by  a  spanwise  bending  mode  whose  deformations  are  moderately  in-phase  with 
the  flapping  motion.  As  the  wing  travels  downward  through  the  midstroke,  the  rigid  body 
velocity  is  augmented  with  a  downward  deformational  velocity,  increasing  the  local  lift  vector. 
The  angle  of  attack  tilts  this  vector  forward,  resulting  in  increased  thrust.  The  phase  lag  between 
the  deformation  and  rigid  flapping  is  large  enough,  however,  to  ensure  that  a  portion  of  the  wing 
is  in  motion  at  all  times,  preventing  the  development  of  drag.  Stanford  and  Beran  note  a  similar 
phenomenon  with  regards  to  thrust-optimal  plunging  shells  . 


Figure  29.  Out-of-plane  wing  deformation  throughout  the  flapping  cycle  of  the  thrust- 

optimal  stiff  topology  design. 


The  power-optimal  design  requires  less  than  half  of  the  peak  power  than  that  of  the  lift-  and 
thrust-optimal  designs,  but  suffers  from  potentially  unacceptable  drops  in  lift  and  thrust,  with  the 
former  particularly  severe.  The  aeroelastic  deformations  that  cause  these  shifts  in  aerodynamic 
loading  are  shown  in  Figure  30.  A  large  adaptive  washout  of  the  trailing  edge  occurs  as  the  wing 
travels  downward  through  the  midstroke,  which  is  a  direct  result  of  the  minimal  support  provided 
by  the  single  diagonal  batten  running  through  the  wing.  The  chordwise  wing  deformation 
attempts  to  align  each  wing  section  with  the  flow  so  as  to  reduce  the  aerodynamic  resistance  and, 
consequently,  the  power  draw.  Adaptive  wash-in  occurs  as  the  wing  travels  upward  through  the 
midstroke,  though  the  trailing  edge  deformation  is  less. 
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Figure  30.  Out-of-plane  wing  deformation  throughout  the  flapping  cycle  of  the  power- 

optimal  stiff  topology  design. 


Though  unsteady  effects  can  lead  to  phase  shifting,  peak  inertial  loads  are  expected  at  stroke 
reversal,  where  rigid  body  wing  acceleration  is  largest.  Peak  aerodynamic  loads  occur  through 
the  midstroke,  where  the  wing  velocity  is  largest.  The  power-optimal  wing  exploits  the  latter,  as 
the  small  deformation  at  stroke  reversal  and  the  sparse  skeletal  topology  both  indicate  small 
inertial  loads.  The  thrust-optimal  design,  however,  appears  to  have  an  equal  presence  of  both 
inertial  and  aerodynamic  loads,  due  to  the  relative  magnitudes  of  deformation  at  the  midstroke 
and  at  stroke  reversal,  as  well  as  the  dense  batten  topology.  The  relative  importance  of 
aerodynamic  and  inertial  loads  in  biological  flight  has  long  been  debated^^.  The  current  results 
would  indicate  that  the  system  might  exploit  either  to  achieve  the  desired  aerodynamic  metrics. 

The  aeroelastic  design  for  the  lift-optimal  configuration  corresponds  to  large  force 
production,  with  the  highest  lift,  thrust,  and  power  history  peaks.  The  deformations  over  a 
flapping  cycle  are  shown  in  Figure  31.  A  large  triangular  portion  of  the  wing  is  composed  of  a 
free  membrane  surface,  constrained  along  most  of  its  border  by  stiff  battens.  This  membrane 
inflates  during  the  downstroke,  increasing  the  camber  and  the  resulting  lift  force.  Similar 
aeroelastic  topology  optimization  strategies  have  been  noted  concerning  the  use  of  perimeter- 
reinforced  membranes  for  lift  production*’®.  As  with  the  thrust-optimal  design,  inertial  loads  don’t 
have  a  significant  role  in  the  optimal  fluid-structure  interaction  of  this  design. 

This  application  combines  a  cellular  division  algorithm  and  a  genetic  optimization  algorithm 
with  an  efficient  aeroelastic  model  to  study  the  effect  of  stiffener  arrangements  on  the 
performance  of  a  flapping-wing  vehicle  in  forward  flight.  Examination  of  the  topologies 
obtained  along  the  Pareto  front  and  their  associated  aeroelastic  behavior  provided  insights  into 
aeroelastic  mechanisms  that  can  be  exploited  through  careful  selection  of  stiffener  layout  to 
improve  performance.  Spanwise  bending,  adaptive  washout,  and  adaptive  cambering  can  be  used 
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to  select  a  desirable  tradeoff  between  thrust  and  lift  production  and  power  consumption.  The 
relative  importance  of  inertial  and  aerodynamic  loads  was  also  seen  to  be  dependent  on  the 
design  objectives. 


Figure  31.  Out-of-plane  deformation  throughout  the  flapping  cycle  of  the  lift-optimal  stiff 

topology  design. 
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CONCLUSIONS 


The  “Physics-based  Design  of  Micro  Air  Vehicles”  laboratory  task  has  produced  several 
unique  studies  into  the  key  system  parameters  and  underlying  physical  interactions  that  drive 
efficient  performance  of  flapping-wing  MAVs.  These  studies  consistently  support  the  conclusion 
that  aeroelastic  and  multidisciplinary  effects  are  indispensible  to  the  performance  of  this  class  of 
aircraft.  Furthermore,  we  have  demonstrated  the  utility  of  optimization  methods  as  a  means  for 
exploring  large  parameter  spaces  to  discover  performance-enhancing  physical  interactions.  Once 
the  optimization  process  identifies  interesting  parameter  combinations,  more  costly,  high  fidelity 
methods  can  be  used  to  verify  the  discovered  phenomenon  and  provide  details  missed  by  the 
lower-fidelity  design  methods. 

The  optimization  strategy  must  be  tailored  to  meet  the  unique  needs  of  each  application.  The 
physical  modeling  should  provide  sufficient  fidelity  to  capture  the  essential  physical  phenomena, 
yet  fast  enough  to  permit  evaluation  at  many  points  in  the  parameter  space  of  interest.  For  this 
reason,  we  rely  on  a  stable  stocked  with  methods  from  a  variety  of  disciplines,  suitable  for  a 
range  of  configurations  and  fidelity.  Methods  can  be  drawn  from  this  stable  and  harnessed,  as 
appropriate,  for  each  application.  Throughout,  methods  that  are  tailored  to  exploit  features  of  the 
response,  such  as  time  periodicity,  or  that  provide  analytic  sensitivities  are  favored  for  their 
speed  and  accuracy. 

We  have  applied  the  optimization  strategy  to  problems  involving  small-scale,  flapping 
aircraft  operating  in  either  hover  or  forward  flight.  These  have  focused  on  key  interactions 
between  the  unsteady  aerodynamics  and  underlying  structure  and  topology.  The  first  application 
sought  out  beneficial  combinations  of  shape,  structure,  and  kinematics  for  a  flapping  wing  in 
hover.  This  work  demonstrates  that  peak  power  can  be  reduced  by  carefully  controlling  the 
relative  contributions  to  total  power  from  aerodynamics,  inertia,  and  structural  deformation.  The 
second  application  produced  insight  into  appropriate  topologies  for  compliant  mechanisms  that 
might  be  used  to  actuate  wing  flapping.  Results  indicate  that  the  flapping  actuator  cannot  be 
isolated  from  the  aeroelastic  system  describing  the  wing;  they  should  be  designed  together.  The 
third  application  studied  the  impact  of  stiffener  arrangement  in  a  membrane  wing  on  various 
performance  indicators.  This  study  reveals  that  aeroelastic  tailoring  can  be  used  to  vary  the 
relative  importance  of  aerodynamic  and  inertial  loads  to  achieve  a  desired  tradeoff  of  lift  and 
thrust  production  and  power  consumption. 

Our  goal  for  follow-on  work  is  to  build  a  detailed  and  quantitative  understanding  of  the 
influence  of  system  configuration  and  flight  atmosphere  on  the  controllability  and  agility  of 
small,  uninhabited,  fixed-  and  flapping-wing  aircraft.  This  will  involve  a  combination  of 
experimental  and  computational  research.  The  experimental  work  will  focus  on  the 
characterization  of  gust  and  turbulence  and  their  effects  on  small  aircraft  operating  in  close  and 
urban  environments.  The  computational  research  will  build  on  our  strategy  of  using 
multidisciplinary  optimization  as  a  means  of  exploring  large  parameter  spaces  characterized  by 
complex,  nonlinear  interactions.  These  computations  will  involve  coupled  aerodynamics, 
structure,  controls,  and  kinematics. 
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